<Name>
<Class>
<Date>
from matplotlib import pyplot as plt
from scipy.io import wavfile
from scipy import signal
import numpy as np
import IPython
from scipy import fftpack
import math
import imageio
plt.rcParams["figure.dpi"] = 300 # Fix plot quality.
plt.rcParams["figure.figsize"] = (12,3) # Change plot size / aspect (you may adjust this).
class SoundWave(object):
"""A class for working with digital audio signals."""
# Problem 1.1
def __init__(self, rate, samples):
"""Set the SoundWave class attributes.
Parameters:
rate (int): The sample rate of the sound.
samples ((n,) ndarray): NumPy array of samples.
"""
self.rate = rate
self.samples = samples
# Problems 1.1 and 1.7
def plot(self, dft = False):
"""Plot the graph of the sound wave (time versus amplitude)."""
#plot regardless of input
ax1 = plt.subplot(211)
domain = np.linspace(0, len(self.samples)/self.rate, len(self.samples))
ax1.plot(domain, self.samples)
ax1.set_ylim(-32768, 32767)
plt.xlabel("Seconds")
plt.ylabel("Magnitude")
if not dft:
plt.show()
#Only plot if dft is True
if dft:
#calculate
c_k = abs(fftpack.fft(self.samples))
n = len(c_k)
k_indices = np.arange(n)
frequency = k_indices * self.rate/n
#plot
ax2 = plt.subplot(212)
ax2.plot(frequency, c_k)
ax2.set_xlim((0,self.rate/2))
plt.show()
# Problem 1.2
def export(self, filename, force=False):
"""Generate a wav file from the sample rate and samples.
If the array of samples is not of type np.int16, scale it before exporting.
Parameters:
filename (str): The name of the wav file to export the sound to.
"""
if force or self.samples.dtype != np.int16:
#Scaling
scaled_samples = np.int16((self.samples * 32767.0)/max(abs(self.samples)))
wavfile.write(filename, self.rate, scaled_samples)
else:
wavfile.write(filename, self.rate, self.samples)
# Problem 1.4
def __add__(self, other):
"""Combine the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to add
to the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the combined samples.
Raises:
ValueError: if the two sample arrays are not the same length.
"""
if len(self.samples) != len(other.samples):
raise ValueError("sample lengths not the same size")
else:
#add matricies elementwise
new_samples = self.samples + other.samples
return SoundWave(self.rate, new_samples)
# Problem 1.4
def __rshift__(self, other):
"""Concatentate the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to concatenate
to the samples contained in this object.
Raises:
ValueError: if the two sample rates are not equal.
"""
if self.rate != other.rate:
raise ValueError("rates are not equal")
else:
#concatenate matricies
new_samples = np.concatenate((self.samples,other.samples))
return SoundWave(self.rate, new_samples)
# Problem 2.1
def __mul__(self, other):
"""Convolve the samples from two SoundWave objects using circular convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
if self.rate != other.rate:
raise ValueError("Rates not equal")
else:
if len(self.samples) < len(other.samples):
self.samples = np.pad(self.samples,(0,len(other.samples) - len(self.samples)),'constant')
elif len(self.samples) > len(other.samples):
other.samples = np.pad(other.samples,(0,len(self.samples) - len(other.samples)),'constant')
#convolve arrays
convolution = fftpack.ifft(fftpack.fft(self.samples) * fftpack.fft(other.samples))
return SoundWave(self.rate, convolution)
# Problem 2.2
def __pow__(self, other):
"""Convolve the samples from two SoundWave objects using linear convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
if self.rate != other.rate:
raise ValueError("Rates not equal")
else:
n = len(self.samples)
m = len(other.samples)
length = n + m - 1
a = np.log2(length)
a = math.ceil(a)
self_sam = np.pad(self.samples,(0,2**a - n),'constant')
other_sam = np.pad(other.samples,(0,2**a - m),'constant')
#convolve arrays
convolution = fftpack.ifft(fftpack.fft(self_sam) * fftpack.fft(other_sam))
#take first n+m-1 elements
convolution = convolution[:length]
return SoundWave(self.rate, convolution)
# Problem 2.4
def clean(self, low_freq, high_freq):
"""Remove a range of frequencies from the samples using the DFT.
Parameters:
low_freq (float): Lower bound of the frequency range to zero out.
high_freq (float): Higher boound of the frequency range to zero out.
"""
n = len(self.samples)
k_low = math.floor(low_freq * n / self.rate)
k_high = math.ceil(high_freq * n / self.rate)
dft = fftpack.fft(self.samples)
dft[k_low:k_high] = 0
dft[n - k_high: n - k_low] = 0
self.samples = fftpack.ifft(dft).real
SoundWave.__init__().SoundWave.plot().scipy.io.wavefile.read() and the SoundWave class to plot tada.wav.rate, samples = wavfile.read("tada.wav") #pull info from file
a = SoundWave(rate,samples)
a.plot()
SoundWave.export().export() method to create two new files containing the same sound as tada.wav: one without scaling, and one with scaling (use force=True).IPython.display.Audio() to embed the original and new versions of tada.wav in the notebook.IPython.display.Audio(filename = 'tada.wav') #read file for comparison
rate, samples = wavfile.read("tada.wav")
tada_ = SoundWave(rate,samples)
tada_.export('tada_not_scaled.wav') #export file, should be same as previous
IPython.display.Audio(filename = 'tada_not_scaled.wav')
tada_.export('tada_scaled.wav', True) #export scaled file, should be louder than previous
IPython.display.Audio(filename = 'tada_scaled.wav')
generate_note().generate_note() to create an A tone that lasts for two seconds. Embed it in the notebook.def generate_note(frequency, duration):
"""Generate an instance of the SoundWave class corresponding to
the desired soundwave. Uses sample rate of 44100 Hz.
Parameters:
frequency (float): The frequency of the desired sound.
duration (float): The length of the desired sound in seconds.
Returns:
sound (SoundWave): An instance of the SoundWave class.
"""
Rate = 44100 #code from textbook
num_samples = duration * Rate
domain = np.linspace(0,duration, num_samples)
samples = np.sin(2*np.pi*domain*frequency)
return SoundWave(Rate, samples)
a = generate_note(440, 2) #A note
IPython.display.Audio(rate=a.rate, data=a.samples)
SoundWave.__add__().SoundWave.__rshift__().a = generate_note(440, 3)
c = generate_note(523.25, 3)
e = generate_note(659.25,3)
cord = a+c+e #notes added elementwise
IPython.display.Audio(rate=cord.rate, data=cord.samples)
a = generate_note(440, 1)
c = generate_note(523.25, 1)
e = generate_note(659.25,1)
arpeggio = a>>c>>e #notes concatenated
IPython.display.Audio(rate=arpeggio.rate, data=arpeggio.samples)
simple_dft() with the formula for $c_k$ given below.np.allclose() to check that simple_dft() and scipy.fftpack.fft() give the same result (after scaling).def simple_dft(samples):
"""Compute the DFT of an array of samples.
Parameters:
samples ((n,) ndarray): an array of samples.
Returns:
((n,) ndarray): The DFT of the given array.
"""
n = len(samples) #code from textbook
m = np.arange(n).reshape(n,1)
W = np.exp((-2j*np.pi/n)*m@m.T)
return W@samples/n
"""
n = 8192 #power of 2
samples = np.random.random(n)
simple = simple_dft(samples)
sci = fftpack.fft(samples)
#compare textbook dft with scipy dft
np.allclose(n*simple,sci)
"""
simple_fft().simple_dft(), simple_fft(), and scipy.fftpack.fft().
Print the runtimes of each computation.np.allclose() to check that simple_fft() and scipy.fftpack.fft() give the same result (after scaling).def simple_fft(samples, threshold=1):
"""Compute the DFT using the FFT algorithm.
Parameters:
samples ((n,) ndarray): an array of samples.
threshold (int): when a subarray of samples has fewer
elements than this integer, use simple_dft() to
compute the DFT of that subarray.
Returns:
((n,) ndarray): The DFT of the given array.
"""
#code from textbook
n = len(samples)
if n <= threshold:
return simple_dft(samples)
else:
even = simple_fft(samples[::2])
odd = simple_fft(samples[1::2])
w = np.exp(-2j*np.pi/n * np.arange(n))
m = n//2
first_sum = even + w[:m] * odd
second_sum = even + w[m:] * odd
return 0.5 *np.concatenate([first_sum, second_sum])
'''
n = 8192
samples = np.random.random(n)
#time computation times
%time n*simple_dft(samples)
%time n*simple_fft(samples)
%time fftpack.fft(samples)
'''
SoundWave.plot() so that it accepts a boolean. When the boolean is True, take the DFT of the stored samples and plot (in a new subplot) the frequencies present on the $x$-axis and the magnituds of those frequences on the $y$-axis. Only the display the first half of the plot, and adjust the $x$-axis so that it correctly shows the frequencies in Hertz.a.plot(True)
cord.plot(True)
Use the DFT to determine the individual notes that are present in mystery_chord.wav.
rate, samples = wavfile.read("mystery_chord.wav")
mystery = SoundWave(rate,samples)
IPython.display.Audio(filename = 'mystery_chord.wav')
mystery.plot(True)
# Find the Magnitudes
n = len(mystery.samples)
k_indices = np.arange(n)
freq = (mystery.rate * k_indices) / n
magnitude = np.abs(fftpack.fft(mystery.samples))
# Find the frequencies
sort = np.argsort(magnitude[:n//2])[::-1]
note_1 = freq[sort[0]]
note_2 = freq[sort[1]]
note_3 = freq[sort[2]]
note_4 = freq[sort[3]]
print(note_1, 'A')
print(note_2, 'G')
print(note_3, 'C')
print(note_4, 'D')
The notes are...
SoundWave.__mul__() for circular convolution.tada.wav.tada.wav and the white noise. Embed the result in the notebook.rate_wn = 22050 # Create 2 seconds of white noise at a given rate.
white_noise = np.random.randint(-32767, 32767, rate_wn*2, dtype=np.int16)
noise_sound = SoundWave(rate_wn, white_noise)
IPython.display.Audio(rate=noise_sound.rate, data=noise_sound.samples)
#convolution
convo = tada_*noise_sound
#append convolution to itself
new_convo = convo >> convo
IPython.display.Audio(rate=new_convo.rate, data=new_convo.samples)
SoundWave.__pow__() for linear convolution.CGC.wav and GCG.wav using SoundWave.__pow__() and scipy.signal.fftconvolve().SoundWave.__pow__() and scipy.signal.fftconvolve() sound the same.rate, samples = wavfile.read("CGC.wav")
cgc_sound = SoundWave(rate,samples)
IPython.display.Audio(rate=cgc_sound.rate, data=cgc_sound.samples)
rate, samples = wavfile.read("GCG.wav")
gcg_sound = SoundWave(rate,samples)
IPython.display.Audio(rate=gcg_sound.rate, data=gcg_sound.samples)
new_tone = cgc_sound ** gcg_sound
IPython.display.Audio(rate=new_tone.rate, data=new_tone.samples)
scipy_tone = signal.fftconvolve(cgc_sound, gcg_sound)
#IPython.display.Audio(rate=new_tone.rate, data=new_tone.samples)
%time cgc_sound ** gcg_sound
%time signal.fftconvolve(cgc_sound, gcg_sound)
Use SoundWave.__pow__() or scipy.signal.fftconvolve() to compute the linear convolution of chopin.wav and balloon.wav.
Embed the two original sounds and their convolution in the notebook.
rate, samples = wavfile.read("chopin.wav")
choping = SoundWave(rate,samples)
IPython.display.Audio(rate=choping.rate, data=choping.samples)
rate, samples = wavfile.read("balloon.wav")
balloon = SoundWave(rate,samples)
IPython.display.Audio(rate=balloon.rate, data=balloon.samples)
choping_balloon = choping ** balloon
IPython.display.Audio(rate=choping_balloon.rate, data=choping_balloon.samples)
SoundWave.clean().noisy1.wav by filtering out frequencies from $1250$-$2600$ Hz. Embed the original and the cleaned versions in the notebook.noisy2.wav. Embed the original and the cleaned versions in the notebook.rate, samples = wavfile.read("noisy1.wav")
noisy1 = SoundWave(rate,samples)
noisy1.clean(1250, 2600)
IPython.display.Audio(rate=noisy1.rate, data=noisy1.samples)
rate, samples = wavfile.read("noisy2.wav")
noisy2 = SoundWave(rate,samples)
noisy2.clean(1250, 4500)
IPython.display.Audio(rate=noisy2.rate, data=noisy2.samples)
vuvuzela.wav by filtering bad frequencies out of the left and right channels individually.rate, samples = wavfile.read("vuvuzela.wav")
vuv1 = SoundWave(rate,samples[:,0])
vuv2 = SoundWave(rate,samples[:,1])
vuv1.clean(200,500)
vuv2.clean(200,500)
vuv = SoundWave(rate, np.vstack((vuv1.samples, vuv2.samples)))
IPython.display.Audio(rate=vuv.rate, data=vuv.samples)
license_plate.png so that the year printed on the sticker in the bottom right corner of the plate is legible.image = imageio.imread("license_plate.png")
plt.imshow(image, cmap = 'gray')
plt.show()
im_dft = fftpack.fft2(image)
plt.imshow(np.log(np.abs(im_dft)), cmap="gray")
plt.show()
The year on the sticker is...2013
mean = np.mean(im_dft)
im_dft[25:45, 90:110] = 350
im_dft[178:184, 328:336] = 330
im_dft[145:155, 225:240] = 320
im_dft[150:250, 300:450] = 350
plt.imshow(np.log(np.abs(im_dft)), cmap="gray")
new_samples = fftpack.ifft2(im_dft).real
plt.imshow(new_samples, cmap = 'gray')
plt.show()